Skip to content

Conversation

@hyeoksu-lee
Copy link
Contributor

@hyeoksu-lee hyeoksu-lee commented Jun 11, 2025

User description

Description

This PR includes an update on initialization method for turbulent mixing layer. Also, Matlab scripts for computing turbulence statistics are added to examples/3D_turb_mixing/turbulence_stat as well as validation results.

Type of change

Please delete options that are not relevant.

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Something else

Scope

  • This PR comprises a set of related changes with a common goal

If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.

How Has This Been Tested?

Ran examples/3D_turb_mixing on Carpenter

Streamwise velocity field at final time
vel1

Vorticity magnitude at final time
omega_mag

Evolution of Reynolds stress
tstep_0

Evolution of TKE budget
tstep_0

TKE budget averaged over self-similar period
avg_self_similar

Momentum thickness growth rate (see run_turbulence.m for references in legend)
growth_rate

Checklist

  • I have added comments for the new code
  • I added Doxygen docstrings to the new code
  • I have made corresponding changes to the documentation (docs/)
  • I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
  • I have added example cases in examples/ that demonstrate my new feature performing as expected.
    They run to completion and demonstrate "interesting physics"
  • I ran ./mfc.sh format before committing my code
  • New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
  • This PR does not introduce any repeated code (it follows the DRY principle)
  • I cannot think of a way to condense this code and reduce any introduced additional line count

PR Type

Enhancement, Documentation, Tests


Description

• Replace linear stability analysis with synthetic turbulence generation method for mixing layer initialization using Guo et al. (2023) approach
• Add comprehensive MATLAB scripts for turbulence statistics computation including Reynolds stress, TKE budget analysis, and validation against reference data
• Update 3D turbulent mixing layer example with higher Reynolds number (320), increased domain size and grid resolution
• Add new configuration parameters mixlayer_perturb_nk and mixlayer_perturb_k0 for synthetic turbulence control
• Remove deprecated mixlayer_domain parameter and associated complex validation requirements
• Implement pseudo-random number generator with supporting constants and helper functions
• Add documentation for new mixing layer perturbation method and turbulence analysis tools
• Update test cases and parameter dictionaries to support new synthetic turbulence approach


Changes walkthrough 📝

Relevant files
Enhancement
6 files
m_perturbation.fpp
Replace linear stability analysis with synthetic turbulence generation

src/pre_process/m_perturbation.fpp

• Removed complex linear stability analysis code for mixing layer
instability waves
• Replaced s_superposition_instability_wave with new
s_perturb_mixlayer function
• Added spectrum-based synthetic
turbulence generation method using Guo et al. (2023) approach

Implemented pseudo-random number generator and helper functions for
perturbation generation

+146/-527
m_initial_condition.fpp
Update function call for mixing layer perturbation             

src/pre_process/m_initial_condition.fpp

• Changed function call from s_superposition_instability_wave to
s_perturb_mixlayer

+1/-1     
m_constants.fpp
Add constants for pseudo-random number generator                 

src/common/m_constants.fpp

• Added constants for pseudo-random number generator implementation

Defined modulus, multiplier, increment, amplifier, and decimal_trim
parameters

+7/-0     
case.py
Update 3D turbulent mixing layer example configuration     

examples/3D_turb_mixing/case.py

• Updated Reynolds number from 50 to 320
• Increased domain size and
grid resolution significantly
• Added grid stretching in y-direction
with new parameters
• Modified time stepping and output settings

+21/-20 
run_turbulence.m
Add MATLAB script for turbulence statistics analysis         

examples/3D_turb_mixing/turbulence_stat/run_turbulence.m

• Added comprehensive MATLAB script for turbulence statistics
computation
• Implements Reynolds stress and TKE budget calculations

Includes plotting functions and comparison with reference data

Provides mixing layer thickness and growth rate analysis

+601/-0 
average_tke_over_self_similar.m
Add TKE budget averaging and validation script                     

examples/3D_turb_mixing/turbulence_stat/average_tke_over_self_similar.m

• Added MATLAB script to compute time-averaged TKE budget over
self-similar period
• Implements interpolation and normalization of
TKE transport, production, and dissipation terms
• Includes plotting
functionality with comparison to reference data from multiple studies

• Contains helper functions for derivative computation and TKE budget
visualization

+118/-0 
Configuration changes
7 files
m_mpi_proxy.fpp
Update MPI broadcasts for new mixing layer parameters       

src/pre_process/m_mpi_proxy.fpp

• Added MPI broadcast for new parameter mixlayer_perturb_nk
• Added
MPI broadcast for new parameter mixlayer_perturb_k0
• Removed MPI
broadcast for deprecated parameter mixlayer_domain

+4/-4     
m_global_parameters.fpp
Add new mixing layer perturbation parameters                         

src/pre_process/m_global_parameters.fpp

• Added new parameters mixlayer_perturb_nk and mixlayer_perturb_k0 for
synthetic turbulence
• Removed deprecated parameter mixlayer_domain

Set default values for new parameters (nk=100, k0=0.4446)

+6/-2     
m_checker.fpp
Simplify mixing layer perturbation validation                       

src/pre_process/m_checker.fpp

• Simplified validation for mixlayer_perturb to only require p > 0

Removed complex validation requirements for model equations, fluids,
and boundary conditions

+1/-5     
m_start_up.fpp
Update namelist for new mixing layer parameters                   

src/pre_process/m_start_up.fpp

• Updated namelist to include new parameters mixlayer_perturb_nk and
mixlayer_perturb_k0
• Removed deprecated parameter mixlayer_domain
from namelist

+2/-2     
case_dicts.py
Update parameter dictionary for new mixing layer options 

toolchain/mfc/run/case_dicts.py

• Added new parameter types for mixlayer_perturb_nk and
mixlayer_perturb_k0
• Removed deprecated parameter mixlayer_domain

+2/-0     
set_user_inputs.m
Add user input configuration for turbulence statistics     

examples/3D_turb_mixing/turbulence_stat/set_user_inputs.m

• Added MATLAB function to configure simulation parameters for
turbulent mixing layer analysis
• Defines fluid properties, Reynolds
number, grid configuration, and time stepping parameters
• Sets up
domain dimensions and stretched grid parameters for y-direction

Configures system indices for multi-fluid and multi-dimensional
simulations

+57/-0   
.typos.toml
Allow TKE acronym in typo checker                                               

.typos.toml

• Added TKE to the list of accepted words in typo checking
configuration

+1/-0     
Tests
1 files
cases.py
Update test case for new mixing layer perturbation method

toolchain/mfc/test/cases.py

• Replaced complex alter_instability_wave function with simplified
alter_mixlayer_perturb
• Updated test case configuration for 3D mixing
layer with new parameters
• Simplified test case setup removing
bubble-related configurations

+33/-53 
Miscellaneous
1 files
golden-metadata.txt
Update test metadata for new build environment                     

tests/D13FDB23/golden-metadata.txt

• Updated test metadata with new build configuration and environment

Changed from Delta cluster to M4Pro local machine
• Updated compiler
versions and CMake configuration

+44/-68 
Documentation
3 files
case.md
Update documentation for new mixing layer perturbation method

docs/documentation/case.md

• Updated documentation for mixing layer parameters
• Removed
mixlayer_domain parameter description
• Added description of new
synthetic turbulence generation method
• Updated requirements for
mixlayer_perturb option

+3/-6     
README.md
Add README for turbulence statistics analysis tools           

examples/3D_turb_mixing/turbulence_stat/README.md

• Added README file describing turbulence statistics analysis tools

Provides references to relevant literature
• Documents file structure
and usage instructions

+23/-0   
references.md
Add turbulence-related academic references                             

docs/documentation/references.md

• Added reference to Guo et al. (2023) for turbulence generation
methods
• Added reference to Michalke (1964) for hyperbolic-tangent
velocity profile instability

+4/-0     
Additional files
9 files
case.py +0/-2     
reference.mat [link]   
golden-metadata.txt +0/-174 
golden.txt +0/-16   
golden-metadata.txt +0/-174 
golden.txt +0/-10   
golden-metadata.txt +0/-174 
golden.txt +0/-18   
golden.txt +12/-12 

Need help?
  • Type /help how to ... in the comments thread for any questions about Qodo Merge usage.
  • Check out the documentation for more information.
  • @hyeoksu-lee hyeoksu-lee requested review from a team and sbryngelson as code owners June 11, 2025 19:32
    @codecov
    Copy link

    codecov bot commented Jun 11, 2025

    Codecov Report

    Attention: Patch coverage is 76.05634% with 17 lines in your changes missing coverage. Please review.

    Project coverage is 44.00%. Comparing base (4864d36) to head (c937bd9).
    Report is 1 commits behind head on master.

    Files with missing lines Patch % Lines
    src/pre_process/m_perturbation.fpp 75.00% 16 Missing and 1 partial ⚠️
    Additional details and impacted files
    @@            Coverage Diff             @@
    ##           master     #879      +/-   ##
    ==========================================
    - Coverage   45.98%   44.00%   -1.98%     
    ==========================================
      Files          68       68              
      Lines       18629    18433     -196     
      Branches     2239     2228      -11     
    ==========================================
    - Hits         8566     8112     -454     
    - Misses       8711     9017     +306     
    + Partials     1352     1304      -48     

    ☔ View full report in Codecov by Sentry.
    📢 Have feedback on the report? Share it here.

    🚀 New features to boost your workflow:
    • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

    @hyeoksu-lee
    Copy link
    Contributor Author

    @sbryngelson Only frontier test fails. I guess something might be wrong with the pseudo-random number generator, but I am still trying to figure out what causes the failure. If you have any idea on what could be possibly wrong, that would be greatly appreciated!

    @sbryngelson
    Copy link
    Member

    @hyeoksu-lee that is indeed weird. if I recall correctly, the Frontier CPU test passes? If that's correct, then it's very strange the GPU case would fail. I notice the failure is for a very small tolerance error. What does that particular test actually look like? Does it meaningfully change the state variables over the ~50 time steps (or whatever it is)? Perhaps it's just not a great test that is very sensitive to compiler optimizations or floating point roundoff.

    integer :: xint, val
    real(wp) :: x

    x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
    Copy link
    Member

    Choose a reason for hiding this comment

    The reason will be displayed to describe this comment to others. Learn more.

    if you think the error in the Cray/Frontier case is due to the RNG, then you can print the random number it generates as a test and see it in the CI output (I think that should work, I believe we print the simulation output of failed tests in the CI output). also: Carpenter is a cray machine, so you can load cray modules there, compile MFC, and see if it passes tests!

    Copy link
    Member

    Choose a reason for hiding this comment

    The reason will be displayed to describe this comment to others. Learn more.

    using Carpenter for this is probably the fastest way to debug. you'll have to change the modules from the ones ./mfc.sh load uses to the Cray ones. let me know if you can't figure this out, I have Carpenter access (i think)

    @hyeoksu-lee
    Copy link
    Contributor Author

    @sbryngelson Thank you for the suggestions! I will try to verify with Carpenter's cray modules. Not sure if I can find all the required modules though. I'll keep you updated

    @sbryngelson
    Copy link
    Member

    @sbryngelson Thank you for the suggestions! I will try to verify with Carpenter's cray modules. Not sure if I can find all the required modules though. I'll keep you updated

    Let me know if that's the case and I'll try to help (put it in Slack).

    @hyeoksu-lee
    Copy link
    Contributor Author

    @sbryngelson Using Carpenter Cray, I found that the PRNG was sensitive to machine precision errors due to the multiplication of a very large number followed by trimming. I’ve updated the PRNG to make it more robust against such precision issues. Test suite passed on both Carpenter GNU and Cray as well as my local macbook, so I believe the Frontier tests should now pass.

    In the meantime, I will have to update the results of example/3D_turb_mixing as its initial condition has been changed and the statistics may change accordingly. I believe the statistics will remain close to the reference, but I want to make sure. This will take a day or two.

    @sbryngelson
    Copy link
    Member

    @hyeoksu-lee ah, good find, very nice! Looks like it passed just fine. FYI the GT (Phoenix) runners are down or acting very strangely right now but I expect them to show fine results. I will mark this PR as a draft PR and you can mark it 'ready to review' once you check your final statistics.

    @sbryngelson sbryngelson marked this pull request as draft June 27, 2025 13:00
    @hyeoksu-lee hyeoksu-lee marked this pull request as ready for review June 28, 2025 01:42
    @hyeoksu-lee
    Copy link
    Contributor Author

    @sbryngelson Statistics looks good. This PR is now ready for review.

    @qodo-merge-pro
    Copy link
    Contributor

    PR Reviewer Guide 🔍

    Here are some key observations to aid the review process:

    ⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪
    🧪 PR contains tests
    🔒 No security concerns identified
    ⚡ Recommended focus areas for review

    Possible Issue

    The new perturbation method uses MPI broadcasts but only generates random numbers on rank 0. This could lead to synchronization issues or deadlocks if MPI communication fails, and the random number generation logic appears complex with potential for numerical instability.

                    if (proc_rank == 0) then
                        call s_generate_random_perturbation(khat, xi, phi, i, y_cc(r))
                    end if
    
    #ifdef MFC_MPI
                    call MPI_BCAST(khat, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
                    call MPI_BCAST(xi, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
                    call MPI_BCAST(phi, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
    #endif
    Code Quality

    The Cholesky decomposition implementation manually checks for small values and sets them to sgm_eps, but sgm_eps is not defined in the visible code. This could cause compilation errors or undefined behavior.

    if (abs(Lmat(1, 1)) < sgm_eps) Lmat(1, 1) = sgm_eps
    Lmat(2, 1) = Rij(2, 1)/Lmat(1, 1)
    Lmat(2, 2) = sqrt(Rij(2, 2) - Lmat(2, 1)**2._wp)
    if (abs(Lmat(2, 2)) < sgm_eps) Lmat(2, 2) = sgm_eps
    Lmat(3, 1) = Rij(3, 1)/Lmat(1, 1)
    Lmat(3, 2) = (Rij(3, 2) - Lmat(3, 1)*Lmat(2, 1))/Lmat(2, 2)
    Lmat(3, 3) = sqrt(Rij(3, 3) - Lmat(3, 1)**2._wp - Lmat(3, 2)**2._wp)
    
    Duplicate Code

    The derivative computation functions f_compute_derivative_1d and f_compute_derivative_3d contain very similar logic for boundary conditions and central differences. This could be refactored to reduce duplication.

    function dfunds = f_compute_derivative_1d(fun,s)
    
        dfunds = zeros(size(fun));    % initialize discrete derivative vector
    
        for j = 1:length(s) % for each index in s
            % if at bottom boundary, use 1-sided derivative
            dfunds(1) = (fun(2) - fun(1)) / (s(2) - s(1));
            % if at top boundary, use 1-sided derivative
            dfunds(end) = (fun(end) - fun(end-1)) / (s(end) - s(end-1));
            % otherwise, if in the interior of the domain, 2-sided derivative
            for i = 2:length(s)-1
                dfunds(i) = (fun(i+1) - fun(i-1)) / (s(i+1) - s(i-1));
            end
        end
    end
    
    % Compute the wall-normal derivative of a discretized function, fun(y)
    function dfunds = f_compute_derivative_3d(fun,s,dir)
    
        dfunds = zeros(size(fun));    % initialize discrete derivative vector
    
        if (dir == 1)
            % Forward difference at start
            dfunds(1,:,:) = (fun(2,:,:) - fun(1,:,:)) / (s(2) - s(1));
            % Backward difference at end
            dfunds(end,:,:) = (fun(end,:,:) - fun(end-1,:,:)) / (s(end) - s(end-1));
            % Central difference for interior
            for i = 2:length(s)-1
                dfunds(i,:,:) = (fun(i+1,:,:) - fun(i-1,:,:)) / (s(i+1) - s(i-1));
            end
        elseif (dir == 2)
            dfunds(:,1,:) = (fun(:,2,:) - fun(:,1,:)) / (s(2) - s(1));
            dfunds(:,end,:) = (fun(:,end,:) - fun(:,end-1,:)) / (s(end) - s(end-1));
            for i = 2:length(s)-1
                dfunds(:,i,:) = (fun(:,i+1,:) - fun(:,i-1,:)) / (s(i+1) - s(i-1));
            end
        elseif (dir == 3)
            dfunds(:,:,1) = (fun(:,:,2) - fun(:,:,1)) / (s(2) - s(1));
            dfunds(:,:,end) = (fun(:,:,end) - fun(:,:,end-1)) / (s(end) - s(end-1));
            for i = 2:length(s)-1
                dfunds(:,:,i) = (fun(:,:,i+1) - fun(:,:,i-1)) / (s(i+1) - s(i-1));
            end
        else
            disp("Wrong dir input");
        end
    end

    Copy link
    Contributor

    Copilot AI left a comment

    Choose a reason for hiding this comment

    The reason will be displayed to describe this comment to others. Learn more.

    Pull Request Overview

    This PR replaces the previous linear-instability initialization of the turbulent mixing layer with a spectrum-based synthetic turbulence generator (Guo et al. 2023), updates all related configuration and broadcast parameters (mixlayer_perturb_nk, mixlayer_perturb_k0), and adds MATLAB tools (and documentation) for turbulence statistics and validation of the 3D mixing layer example.

    • Replace s_superposition_instability_wave with s_perturb_mixlayer and remove legacy linear-stability code.
    • Introduce new parameters mixlayer_perturb_nk and mixlayer_perturb_k0, propagate them through global parameters, MPI proxies, start-up, and namelists.
    • Add MATLAB scripts (run_turbulence.m, average_tke_over_self_similar.m), example README, and reference documentation.

    Reviewed Changes

    Copilot reviewed 22 out of 31 changed files in this pull request and generated 1 comment.

    Show a summary per file
    File Description
    src/pre_process/m_perturbation.fpp Remove linear-stability routines; add s_perturb_mixlayer subroutine and PRNG constants.
    src/pre_process/m_initial_condition.fpp Update initialization to call s_perturb_mixlayer instead.
    src/common/m_constants.fpp Define PRNG constants (modulus, multiplier, increment, etc.).
    src/pre_process/m_global_parameters.fpp Add mixlayer_perturb_nk and mixlayer_perturb_k0 defaults.
    src/pre_process/m_mpi_proxy.fpp Broadcast new mixing-layer parameters; remove deprecated mixlayer_domain.
    src/pre_process/m_start_up.fpp Include new perturbation parameters in namelist.
    src/pre_process/m_checker.fpp Simplify validation rules for mixlayer_perturb.
    toolchain/mfc/run/case_dicts.py Register new parameter types for mixlayer_perturb_nk and mixlayer_perturb_k0.
    toolchain/mfc/test/cases.py Replace legacy instability-wave test with alter_mixlayer_perturb.
    examples/3D_turb_mixing/turbulence_stat/* Add MATLAB scripts, README, and averaging tools for validation.
    examples/3D_turb_mixing/case.py Update 3D example to higher Re, grid stretching, and mixlayer config.
    docs/documentation/case.md Update docs to describe new synthetic turbulence parameters.
    docs/documentation/references.md Add Guo (2023) and Michalke (1964) references.
    .typos.toml Allow “TKE” in typo checker.
    Comments suppressed due to low confidence (3)

    examples/2D_mixing_artificial_Ma/case.py:110

    • Enabling mixlayer_perturb for a 2D case (p == 0) will trigger the PROHIBIT check (p > 0) and cause startup to fail. Please disable or guard this setting in 2D examples.
                "mixlayer_vel_coef": 1.0,
    

    src/pre_process/m_checker.fpp:191

    • The new PROHIBIT rule only checks p > 0 but does not enforce mixlayer_vel_profile = .true.. To avoid misconfiguration, consider adding validation that the velocity profile flag is enabled when perturbations are requested.
            @:PROHIBIT(mixlayer_perturb .and. p == 0, "mixlayer_perturb requires p > 0")
    

    toolchain/mfc/test/cases.py:634

    • The new test generator only covers 3D (len(dimInfo[0]) == 3) cases. If 2D support is intended, add a branch to test mixlayer perturbations in 2D.
        def alter_mixlayer_perturb(dimInfo):
    

    @hyeoksu-lee
    Copy link
    Contributor Author

    hyeoksu-lee commented Jun 30, 2025

    @sbryngelson I guess this is now ready for merging

    @sbryngelson sbryngelson merged commit 5b8e2eb into MFlowCode:master Jul 1, 2025
    40 of 41 checks passed
    prathi-wind pushed a commit to prathi-wind/MFC-prathi that referenced this pull request Jul 13, 2025
    …esults (MFlowCode#879)
    
    Co-authored-by: Copilot <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Spencer Bryngelson <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    Co-authored-by: Hyeoksu Lee <[email protected]>
    @hyeoksu-lee hyeoksu-lee deleted the mixlayer branch July 16, 2025 22:11
    @hyeoksu-lee hyeoksu-lee mentioned this pull request Jul 17, 2025
    4 tasks
    Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

    Development

    Successfully merging this pull request may close these issues.

    2 participants